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Background 
A direction-finding system consists of both a net control station 


and a set of direction-finding sites at widely dispersed geographic 
locaticns. When instructed by the net control station, each site tunes 
its receiver to a specified frequency and, using extremely directional 
antennas, obtains a bearing on the signal it receives. These bearings 
are then reported to the net control] station where they are processed 
to estimate the true location of the emitter (obtain a "fix"). 

Should all the bearings intersect at a single point, the estimate 
is, ot course, obvious. Such a situation rarely occurs in practice 
since error is induced in the direction of the received signal by both 
the equipment at each site and by the propagation medium between the 
sites and the emitter. On the basis of ampirical evidence, these errors 
appear to be approximately normally distributed. Thus, the reported 
bearings may be thought of as statistical samples taken from normal dis- 
tributions with unknown true means. Various statistical techniques may 
be used to estimate these means, thereby providing an estimate of the 
true location of the emitter. 

A particularly vexing aspect of the problem, however, is that of 
so-called “wild bearings." When a site tunes up on the specified 
frequency, it might not hear the same signal that other sites are 
receiving. This can happen both as a result of propagation conditions 


and of slight differences in the times the various sites listen. In the 
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latter case a site might report a bearing on the wrong emitter. Fram a 


statistical point of view, a wild bearing cannot properly be considered 
to be a sample from a normal distribution with maan bearing on the true 
emitter location. Unfortunately, it is Lnpossible to distinguish 
between bearings taken on the correct emitter from bearings taken on 
another emitter prior to processing. All fix computation algorithms 
must deal with this problem in one way or another. 

A common method is to consider the problem in two phases. In the 
first phase, the “most likely" set of true bearings is identified and 
in the second phase only these bearings are used in making the estimate. 
The difficulty with this approach is that in deciding which bearings are 
wild, it is necessary to establish a preliminary estimate of the true 
location. Then bearings which are not close to this position are con- 
sidered wild and are rejected. If the preliminary estimate is not 
Close to the true location, it is quite likely that good bearings will 
be thought to be wild and wild bearings good. It is also possible that 
enough bearings will be rejected so as to create a "no-fix" situation, 


i.e., not enough bearings remaining to reliably compute a fix. 
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Overview of the FALCONFIX Method 


The method developed in this paper takes a considerably different 
approach to the problem-no bearing is ever rejected. Instead, the 
bearings reported fram each site is considered to be a sample from a 
nomal distribution with probability p (a good bearing) and froma 
uniform distributicn with probability l-p (a wild bearing), where p is 
determined empirically by observing the fraction of good bearings over 
a period of time. Then a likelihood function, the joint probability 
distribution of bearing errors assuming this composite normal/uniform 
distribution, is defined over a spherical earth model. The estimate of 
true enitter locacion is then taken as the location which maximizes the 
likelihood. Because of the nature of each of the bearing error distri- 
butions, this likelihood function may have numerous relative maxima 
(modes), corresponding to points an the earth where intersections of the 
reported bearings are concentrated. 

Non-linear programming and regression techniques are then used to 
find the location and value of each mode of the likelihood function. 
The algorithm generates the three most likely modes and crders them 
according to the probability that they are associated with the true 
emitter location. The net control station operator may then determine 
which position to report, an the basis of both the algorithm's ordering 
as well as other available information. 

Finally, the algorithm generates an elliptical confidence region 
with a specified probability of containing the true position of the 


emitter. 


Chapter 2 
MATHEMATICAL DESCRIPTION 
Distribution of Bearing Error 
The bearing which is reported from a direction-finding site may be 
taken on either the true target or, inadvertently, on a false target. If 


it is taken on a true target, we assume that the angular error (¢) between 
the reported bearing and the true bearing is a random variable with a 
normal distribution having mean zero and variance 0°, On the other hand, 
if the bearing is taken on a false target, no information is conveyed 

by the reported bearing concerning the bearing to the true target. 
Therefore, we assume that the angular error is a random variable with a 
uniform distribution. Further, we assume that the bearing is taken on 
the true target with probability p and on the false target with 
probability (1-p); the value for p may be cotained fron empirical 
evidence. For small values of ¢€, therefore, we assume the distributic¢n 
to be normal ard for iazrger values we assume it uniform. The transition 
points are designated +co, where c is a constant which depends upon both 
pando. Lastly, it is assumed that the probability density function 

of the bearing error is a continuous functir: of €. The resulting dis- 


tribution, which we shall call “canposite nomnal-uniform" is represented 


by 
1 exp(-e2/20°) , -coseSor 
Kv 2110 
h{e)= 
= exp(-c/2), -1Se<-co or co<eSm 
K/2n0 


where € and o are in radians and the normalizing constant K is given by 
K = [26(c)-1] + 2[n-co} fev 
4210 
where $ is the standad normal cumulacive distribution function. 
The graph of a typical composite normal-uniform is depicted in 
Figure 1 below. The fact that the ordinate is constant for values of 
€ greater than co or less than ~co is extremely significant to the 


‘ethod, as shall Leccme apparent in later sections. 


Figure 1. Typical Camposite Nomal-Uniformm PDF 
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To determine the transition point, we note that 


P| 1 hase pee 
p-zf = er eae 


1 
=e [26(c)-1), or 


p= 2¢(c)-1 
[20(c)=1} + 2{n-co} ia ae Cet a 


mC 
In Figure 2 below, this relationship between p and o is plotted for 
various values of c. Then, given values of p and a, c can be obtained 


by interpolating between the curves. 
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Figure 2. Transition Point Curves 


One final approximation for the density function is required to 
make the procedure more efficient computationally. The e* tem is 
difficult to work with using vector methods, so the Maclaurin approxi- 
mation for cos € is used tc derive the relationship 

c*m 2(1 - cos €), 
This approximation is correct te within 0.4% for values of c« up tc 0.2 


radians, a value greater than will occur in practice. 


Thus, h(c) is approximated by 


( : on | (1-cos a] , scoSesoc 
oO 


fle) = 


J exp E (1-cos (co) } , ~m™Se<-co or co<eSn 
Oo 
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The Likelihood Function 


The likelihood function for the emitter location (¢,9) is defined 
to be: 
F($¢,8) = GlEyr€ny000r€,) = £) (e,)£, (€5)+++f, (e,) 
where €; is the angular error between the reported bearing fram station 
i and the true bearing from station i to (>,8). Thus or is a randem 


variable with density function 


1 


Sak we <- < 
exp — {1-cos "| 1 “C,I;SE;Sc,0; 


f.(e€.) = 
a 1 1; 
exp |- l-oos (c,0,)) -1S€.<-¢,0, or 
K./2n0, 0.2 ii 5 iit 
: = ‘ c.0.<e,£n 
1s ir 
Therefore, 
n 
To, 8) -—o3 ih ~ it exp =}, (1-cos €.)1X 
(2m) ii 0. = 
isl icW = 


-1 
bi exp [—7(1-cos (c;9;)) 
: 0. 
itw i 


where W= {i: 


< 
oy |Se, 9, } 


As will be made clear subsequently, our approach to estimating 
emitter location will be to find the relative maxima--which we call 


"“modes"--of this function. 
n 
Since L “aa. K.o. is not a function of the c F(¢,8) is 
one OP ir FU, 
isl 


maximized whenever 


1 1 

F*($,9) = il exp |- —5(l-cos €,) ji exp [- —5(1-cos (c, 0, )) 

f o, 
iew i i¢w i 


is a maximm. F*(¢,9), in turn, is maximized whenever its logarithm 


L(¢, 8) -) cos e,) -> 4 (100s (¢;9,)) is, and this is the 
| icew itw i 
| function to which we shall direct our attention. 

Before proceeding with the development of the method, we shall pause 
to consider exactly what this function L(¢,6}) looks like. In Figures 3 
through 7 are depicted five examples of a set of bearings and the 
resulting function L. 

In the figure below the case in which two lines of bearing intersect 
at a single point is illustrated. Notice that the function has a single 
mode at the point of intersection. (As we shall see later, the FALCONFIX 
algorithm does not consider two bearing crosses as potential position 


extimates. } 
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Figure 3. L(¢,9); Two Bearing Cross 


In Figure 4 below, the function L(¢,6) in which three bearings 


intersect at a point is illustrated. Here the mode is samewhat more 


sharply peaked than in the previous case. 
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Figure 4. L($,8); Three Bearing Case #1 


In Figure 5, the three bearings do not meet at a single point, and, 
in fact, the "cocked hat" is so large compared to the standard 
deviation of the bearing errors, that the function does not have a 
single mode in the center. Instead, there are three individual two 


bearing crosses at each of the intersections. 
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Figure 5. L($,6); Three Rearing Case #2 


In Figure 6 below, a four bearing case is depicted in which the 
bearings do not intersect at a single point. Notice that this figure is 
considerably less regular than the previous figures, and is more typical 


of what one finds in real cases. 
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Figure 6. L($,8); Four Bearing Case 


In Figure 7 the bimodal case is illustrated. In this five bearing 
example it is not likely that all bearings were taken on the same 
emitter. Rather, it is much more plausible that either Bearings 1, 2 
and 3 are good bearings (with 4 and 5 wild) or Bearings 1, 4 and 5 are 
good bearings (with 2 and 3 wild). The two modes in the figure corres- 
pond to these two possibilities. 


-_ ae 


Figure 7. L($,6); Bimodal Case 


Before going on to the methodology of identifying these modes, it 
will first be necessary to express each of the ej in the equation for 
L explicitly as functions of ¢ and 6, the latitude and longitude of 
the position estimate. To accemplish this we will rely heavily on the 


use of vectors as described in the next section. 
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Vector Representation 


A position vector P from the center of the earth to a point ($,9) 
on the surface may be represented as 
P = (cos 8 cos $)i + (sin 8 cos $)j + (sin $)k 
This vector P is depicted in Figure 6 below in a right-hand 
coordinate system in which I points to the prime meridian on the equator 
and k points to the north pole. Thus, east longitude is represented by 
positive values of 6, west longitude by negative values, north latitude 


by positive values of ¢ and south latitude by negative values. 


Figure 8. Position Vector 
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If (¢,8) is the location of the target and (¢, 685) is the location of 
the it? station, then the corresponding position vectors are: 
T = (cos 6 cos ¢)i + (sin 8 cos $)j + (sin ¢)k 
S, = (cos 6, cos 9,)i + (sin 8; cos ¢)5 + (sin ¢,)k 
Let 8, be the bearing reported by the i‘? station. Then in terms of 
the local coordinate system (i*, j*, k*) in which I* points east, }* 
points north and k* points to the local zenith, a vector pointing in 
the direction of Bs may be represented as 
B,* = (sin B.)i* + (cos B.)j* 


which is shown in the figure below. 


Figure 9. Bearing Vector 


To be useful, however, it is necessary to express each bearing 
vector in the same coordinate system as the T and S; vectors. This is 
accamplished by rotating the local frame with transformation R; so 


that it coincides with the geocentric frame shcwn in Figure 6. 
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Specifically, 
2 e = z 5 k 
B, R,B; Ai + B} + Ck where 
A=- {sin ¢, cos 8; cos B, + sin 8; sin B,) 
B = oos 6; sin 8, ~ gin $, sin 8, cos Bs 


C = cos $; cos B. 
Now, T, the S. and the B, are all unit vectors expressed in the 


same geometric frame. 


S.XxT 
Then, the vector is a unit vector normal to the "station- 
IS xt| 
Sx Be 
target plane," the plane containing Si and T. Similarly, aa 
|S.xB. 
5 Wes & 


is a unit vector normal to the "sStation-bearing plane," the plane con- 


taining S; and 3. 


Clearly, the angle between these two normal vectors is equal to 
the dihedral angle between the two planes, which in turn is equal to 
the angle ey between the reported bearing and the true bearing from 
the Aa station to the target. This relationship is shown in Figure 10. 


Fran the familiar dot product relationship we have 


Rw pe RR er ee 
e 


S xT xB, \ 


; 7 - 
IS.xT| 59,1] 

Since S, and B, are orthononnal, it is true that |s,x B,| = 1, which 

means tnat 

(S XT) . (S > B) 


oos €, = 
[Sixt 
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Figure 10. Dihedral Angle 
Employing a little-used identity frem vector calculus on the numerator 
of this expression yields 


But S.°S; = ] since 8; is a unit vector and S,°B; = 0 since S; and By 


are orthogonal. 


Thus (S,xT) + {S XB) = B;°T, which means that 


B.°T 


COS ©. = 
* [sx] 


Now [s.xt| is equal to the sine of the angle between S. and T. 
Since S,°T equals the cosine of the same angle, we can represent 


|sxT] as 


Isat = \fn is -1)2 
which is somewhat easier to deal with computationally. The relationship 


B.°T 
cos ¢. 2 —o 
1 


Vi1- (S;-T) 


will be used in subsequent calculations. Making this suostitution, the 
logarithm of the likelihood function becomes: 


-1 Be ai 
L(¢, 8) >) 5 1-— . mer } (1-cos (c;0,)) 
iow “i Vi-(s,-7)? itw 7d 


In this equation the Cir Oy, By and S. are all constants and T is 
a function of ¢ and 8. We now have L(¢,6) in a form where it is 
possible to find the modes by moving T in such a way as to progressively 


increase the value of L. In the next section we will illustrate exactly 


hew this is done. 


Modes of L(¢, 6) 
L(¢,8) can be written in the fom: 


B,°T 
L(o,8) -) os ek 


iew 0,7 V1- (S;-T) 


where K = a constant 
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Be = bi ji + b; 9) + b. 3k 
S; = $31 + S; J + S; 3k 
T= t,2 + t,J + t5k 


cos 6 cos ¢ i+ sin 6 wos ¢ 3 + singk 
and 'B,! = |S,| = |T| =1 
A necessary condition fc. the point (¢*,6*) to be a mode (relative 
maximum) of L(¢,9) is that both 


au! 
36 | o=9* 


After considerable calculation it can be shown that: 


aL = _ 
$6 |o=o* =0 and = 0 


~b, ,sin8 cos} + b. cos® cosd 
% 5 F. (7) il i2 


30 i BT + 
icW 
(S;°T) (-s, )siné cos + 8; ,c0s8 cosé) 
+], ara, 
1-(S. +) 
pA 
aL j i is ~b, ,cosb a - b, ,siné sind + b, soos Pe 
ad i B;°T 
iew 


(S;-T) (-s, ;cos6 sind - S; 28in8 sind + 8, 0080) 


2 
1- (S; -T) 


— 


B,°T 
where F; (T) = 


2h (5, +7)" 


%4 4 
Setting these partial derivatives equal to zero and then solving 


for 6 and > respectively, we obtain: 


bio (S;+T) Si> 


Py eo aaet 
-1 Jiew : = 1-(s,°7)? 
§ = J(o,8) = tan - 
>F is pts la 
Pe" Bet ig. ony? 
ieW i 
and 
bB: (S,°T)s. 
i3 i i3 
> F, (7) ST + oes 
1 iew 1 ~(S;+T) 


¢ = K(¢,0) = tan 


vy, = b, ;cos8+b; »siné i (S;-T) (s, costs, osine) 
& i B,°T 1- (5, +7)? 

Thus, at the point (~*,9*) these equations will be satisfied. These 
equations are used to establish the following recursion relationships: 
Shel = J (9,8) ) and Sie = K(9, 844) , Kh 220,12, 8). 

The FALCONFIX algorithm for estimating (¢*,9*) proceeds as follows: 
l. Select a starting point ($5280) . (The selection procedure 
is discussed later.) 


2. Tp = T ($5089) : Initialize the target vector 


3. k=0 : Initialize k 


4. 6 = 6) : Specify the stopping rule 


5S. 94) = 70,8) ; Calculate the next longitude 

6. oy = KH 87) : Calculate the next latitude 

7 Try = TU 41-841) + Calculate the next target vector 

8. If tom then 6=6, : Relax, the stopping rule if k is large 
9. Tf |T,)-T,1<6 or k>N then go to 12: stopping rule 

10. k = k+l : Increment k 

ll. Goto 5 : Perform the next iteration 


‘ : i 
12, Cree is an estimate for (*,8*) 


Note: N>M and $7 5) 


The intersections of pairs of bearings are used as the starting 
points for the algorithm. Those pairs which intersect at an angle of 
greater than 20° are considered first—those intersecting at smaller 
angles are deferred for later consideration. The first check in deter- 
mining whether an intersection constitutes a valid starting point is if 
that particular pair is not a subset of a bearing set W from a previously 
generated mode. If it is, the intersection is rejected since the 
iteration procedure would yield the same mode. The second check is to 
determine if at least cone additional bearing is contained in the set W 
at the point of intersection. If not, the intersection is also rejected 
as a starting point since in this case no iteration will take place and 
two-bearing modes are not allowable, except in the special case where 
only two bearings are input. 

The following example illustrates the algorithn and the selection 


of starting points. Suppose the reported bearings are as shown in 
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Figure 11. Example of Mode Determination 


Figure 11. The incersection (1,2) is temporarily set aside because the 
angle is too small. (1,3) satisfies this criterion and also that at 
least one additional bearing, (2), is contained in the bearing set at 
this point. Therefore (1,3) is a valid starting point. The algorithm 
camences and after several iteratic... converges to the point A. Since 
bearings 1, 2 and 3 are in this bearing set, the intersections (1,2), 
(1,3) and (2,3) are removed fram further considerations. (1,4) is the 
next. candidate and since it meets both criteria it is selected as a 
starting point. The algorithm begins there and converges to point B. 
In the same way as before, (1,2), (1,4) and (2,4) are eliminated. 


(3,4) is the only intersection remaining. It satisfies the angular 


criterion, but no other bearing is close eregh to be a member of the 


bearing set at the point of intersection. Therefore (3,4) is eliminated. 
Since no intersections remain as candidates, all modes have been faurd 
and the problem is terminated. 

When the algorithm generates more than a single mode, we order them 
with respect to the probability that they coincide with the true location 
of the emitter. Since each bearing has a probability p of being a true 
bearing and for realistic direction-finding systems, this value is 
greater than 0.5, the mode having the greatest number of forward bearings 
in its bearing set is considered the most likely. In the case where two 
or more modes involve the same number of forward bearings, the algorithm 
uses the degree of dispersion of the bearings as a secondary criterion. 
Specifically, the average squared perperdicular distance from the position 
estimate to each line of bearing is computed and the modes are ordered 


so that the one with the smallest such value is first. 
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Standard Deviation of Bearing Error 
The standard deviation of the bearing error (9) is used in the cal- 


culation of the Likelihood function and its partial derivatives. While 
the particular method used for obtaining 0 is not crucial to the 
FALCONFIX algorithm, for the sake of completeness its computation will 
be discussed. 

The standard deviation is taken to be a function of both the charac- 
teristics of the particular receiving site (antenna, equipment, location, 
exc.) am! of the distance between the site and the emitter location. In 
particular, we assume the multiplicative relationship o = h(d)o where o 
is a constant representing the normalized standard deviation of the 
bearing error from the site and h(d) is an adjustment for distance, 

If s = sin §, where § is the central angle between the target vector 
T and the station vector S, then 


d= 9444ecin ¢(s) , d in nautical miles 


or 


s = sin ed 
3444 


The function h(d) below is motivated by the Ross Range Curve, 


although certain moditications have been made, 


60 , if ds 60 
60-3152. 387 (s-.0174524), if 60<ds 120 
naire -174524/s , if 120<ds 600 
s + .826352 , 1f 6005S ds 5410 
1.826352 , if d> 5410 


A graph of h(d) is provided in Figure 12 below, 
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DISTANCE (NAUTICAL MILES! 


Figure 12. Graph of h(d) 

The extremely large standard deviation for distances less than 100 
nautical miles or so, is essentially an artifice to assure inclusion of 
bearings fram sites which are close to a potential mode; such bearings 
are almost invariably good bearings. For example, if 5 is equal to 1.0 
degrees and if the transition point c (referred to in Figure 2) is set 
equal to 3.0, then any estimate of emitter location within the “keyhole” 
in Figure 13 will include the bearing fran site s. (In this figure the 
origin represents the site and the horizontal axis is the direction of 
the reported bearing). Notice that the bearing is included whenever the 


site is within 60 nautical miles of the estimate. 
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DISTANCE «NAUTICAL MILES: 


Figure 13, Inclusion Region 
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The REFIX Procedures 


One final procedure is used to refine the position estimates 


generated by the previously described algorithm. The need for such a 


termination procedure is two-fold. 


- If one of more of the bearings in the included set is from a 
station which is reasonably close to the position estimate, the 
artificially large standard deviation associated with such 
bearings will tend to diminish its rightful impact on the 
position estimate. 

- If the likelihood function happens to be unusually flat in the 
neighborhood of the mode, it is possible that the stopping rule 
may terminate the algorithm prematurely, resulting in a reported 
position estimate which is not coincident with the point at 


which the relative maximm occurs. 


Two different termination procedures have been developed, but both 


have the following features in common. 


~ The bearing set is "frozen" to include only those bearings 
which were included in the final iteration of the algorithm 
previously described. 

- The standard deviations of bearings fram stations within 600 
nautical miles of the position are modified by changing h(d) to 


h* (d) = 0.75 + a@/2400 where d is the distance in nautical miles. 


REFIX/1. In this procedure, after the bearing set has been frozen 


and the standard deviations of bearings fran close-in stations have been 


modified, the iterative algorithm is employed again. 


REFIX/2. This termination procedure uses techniques employed in 
multiple regression analysis. Most of the computation done here is 
required for the confidence region calculations which are discussed in 
the next section. In this procedure a cartesian coordinate system is 
set up in a plane tangent to the spherical eartn model, such that the 
origin is coincident with the position estimate produced by the iterative 
procedure. The lines of bearing are assumed to be straight lines in this 
plane. 

A single representative line of bearing is shown in Figure 14 with 


(XY) being the true target location. 


Figure 14. Coordinate System for the REFIX/2 Procedure 


Fron elementry algebra, the equation of the solid line in Figure 14 


is (-cos 8;)x + (sin By)y = &;. The equation of the dashed line through 


the true target location (x, +¥,) is (-cos 8.)x + (sin Boy 2 é -e@ 
Then, since (x,-¥,) is on the dashed Line, we have: 


Sin (-cos B.)X, + (sin Boy, ~ @; 


This equation is in the form of a general multiple regression model. 
Both the 8 and the 6. can be calculated since the origin is known, The 
Bs is approximately the bearing from the origin to the station for 
distant stations and approximately the reported bearing for near stations. 
The 8; is approximately v5 sin qd; where e is the angular difference 
between the reported bearing and the becring fram the station to the 
origin, and where qd; is the distance from the station to the origin. 
Also, en the perpendicular distance from the true target location to 
the bearing, is a normal random variable since e, Fe, sin qd. and e5 
(the angular error) is a normal random variable. Furthermore, since the 
€,; are normal and independent with mean zero and variance a.” the e; 
will be normal and independent with mean zero and variance (sin“d,)o,7. 
_ Thus, the problem reduces to the simple multiple regression problem of 
estimating b when 


S=Abt+e 


eee ee Pe 


Further, we have that e is distributed as a multivariate randon 


variable with ) = Qed and diagonal covariance matrix 


sin?d w 2 Qo * 8 *¢ re) 
Vil 22 
re] sin, WS 
mi fed 2 . ry . = a 2 wW 
P= 9% ; . , re) 
re) a Sa | oe sin’dw,? 
with a7 = W554) for each station. 


Using the above notation, the maximm likelihood estimate of 


b= (X56¥,) the true target location, is given by 


bea clattal rls 
Once b is available, we update the target vector to give the refined 
position estimate. 

Relative Advantages. Either termination procedure works quite 
adequately in practice. The REFIX/1 procedure has a slight advantage 
when the position estimate is very close to one station since the 
assumption of parallel movement of bearings which is made in REFIX/2 is 


not really true in this case. On the other hand, REFIX/2 is quicker in 


some cases since on occasion REFIX/1 requires a number of additional 


iterations. 


The Confidence Region 


The calculation of the confidence region follrvs directly from the 


REFIX/2 algorithm. The estimate 5 is distributed multivariate normal 
with mean U* = hey)" and covariance matrix 
ree at stabs atwhay? 
°o 
The minimm area 100 (l-a)® confidence region is given by 


a x 
—z ley) E . FA 2ynee 


s 


NO 


‘ 2 . ; 
where s, is the estimate of 0, based on n data points and Fy ,2,m2 in 
the appropriate F distribution percentage point. In this algorittm, 


the 0,” is assumed known so Fo 5% is used. This value is approximated 


by ~log, (a). 
The ellipse which is generated is given by the direction of the 


major avis and the magnitudes of the semi-major and semi-minor axes. 
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Computer Program Flow Chart 
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Camputer Program Listing 
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